////////////////////////////TCF/////////////////////////////
// TPC Coherent Fit (TCF)                                 //
// Author: Xianguo Lu, lu.xianguo@gmail.com               //
// Reference: Xianguo Lu, "Exploring the performance      //
// limits of the ALICE Time Projection Chamber and        //
// Transition Radiation Detector for measuring identified //
// hadron production at the LHC", dissertation at the     //
// University of Heidelberg (2013).                       //
////////////////////////////////////////////////////////////

#include "MN2D.h"

ClassImp(MN2D);

Int_t MN2D::fkParaResolution = 0;
Int_t MN2D::fNPolint = -999;
Double_t MN2D::fTolerance = -999;

Int_t MN2D::fRegOpt = 0;
Bool_t MN2D::fkFixMean = kFALSE;
Double_t MN2D::fElectronHighPTrigger=-999;
Double_t MN2D::fElectronHighPFrac=-999;
Double_t MN2D::fRegWeight = -999;
Double_t MN2D::fCpIni = -999;

Bool_t MN2D::fkUseTrueYield = kFALSE;

TH2D *MN2D::fHoriginal=0x0;
TH2D *MN2D::fHFitData=0x0;
TH2D *MN2D::fHprefilter=0x0;
TH2D *MN2D::fHrest=0x0;
TH1D *MN2D::fHFitYield=0x0;
TH1D *MN2D::fHTrueYield=0x0;
TH1D *MN2D::fHrestyield=0x0;

Int_t MN2D::fNbin = 0;
Int_t *MN2D::fBinX = 0x0;
Double_t *MN2D::fXCenter = 0x0;
Double_t *MN2D::fFitYieldArray = 0x0;
Double_t *MN2D::fTrueYieldArray = 0x0;

Double_t MN2D::fMinXing = -999;
Double_t MN2D::fMaxXing = -999;

Bool_t MN2D::fkPreID = kFALSE;
Double_t *MN2D::fPrePar = 0x0;
Double_t *MN2D::fPreAll = 0x0;

Int_t MN2D::fNYbins = 0;
Double_t MN2D::fYbw = -999;
Double_t *MN2D::fYCenter = 0x0;
Int_t MN2D::fNdata = 0;

Int_t MN2D::fNfix = 0;
Bool_t *MN2D::fIfixed=0x0;

Double_t *MN2D::fLows=0x0;
Double_t *MN2D::fHighs=0x0;
Double_t *MN2D::fPars=0x0;
Double_t *MN2D::fLowErrs=0x0;
Double_t *MN2D::fHighErrs=0x0;
Double_t *MN2D::fAs=0x0;
Double_t *MN2D::fErrAs=0x0;
Double_t *MN2D::fRatio=0x0;
Double_t *MN2D::fErrRatio=0x0;
Double_t *MN2D::fACov=0x0;

Double_t MN2D::fStatChi2=-999;
Double_t MN2D::fBeta = -999;
Double_t MN2D::fPenalty=-999;

FILE * MN2D::fLog = 0x0;
TStopwatch *MN2D::fStopwatch = 0x0;
TMinuit *MN2D::fMinuit = 0x0;
Int_t MN2D::fkErr = -999;

//==============================================================================
//                                    initialization and destruction
//==============================================================================

void MN2D::Ini(TH2D * hpreID, TH2D * horiginal, const Int_t nb, const Int_t binx[], const Double_t bbpar[], const Double_t binpar[], const Double_t realRes)
{
  const TString pwd = gSystem->pwd();

  TPCHist tpchist(horiginal);

  //================================== fkParaResolution
  fkParaResolution = -999; 
  if(pwd.Contains("RES1000")){
    fkParaResolution = 1000;
  }
  else if(pwd.Contains("RES2000")){
    fkParaResolution = 2000;
  }
  else if(pwd.Contains("RES3000")){
    fkParaResolution = 3000;
  }
  else{
    printf("MN2D::Ini wrong fkParaResolution %s\n", pwd.Data()); exit(1);
  }

  //======================================== fNPolint

  fNPolint = -999;

  if(pwd.Contains("POLINT-999")){
    fNPolint = 10000;
  }
  else if(pwd.Contains("POLINT10")){
    fNPolint = 10;
  }
  else if(pwd.Contains("POLINT14")){
    fNPolint = 14;
  }
  else if(pwd.Contains("POLINT18")){
    fNPolint = 18;
  }
  else if(pwd.Contains("POLINT22")){
    fNPolint = 22;
  }
  else if(pwd.Contains("POLINT2")){
    fNPolint = 2;
  }
  else if(pwd.Contains("POLINT4")){
    fNPolint = 4;
  }
  else if(pwd.Contains("POLINT6")){
    fNPolint = 6;
  }
  else if(pwd.Contains("POLINT8")){
    fNPolint = 8;
  }
  else{
    printf("MN2D::Ini wrong fNPolint! %s\n", pwd.Data());exit(1);
  }
  //========================== fTolerance

  fTolerance = -999;
  if(pwd.Contains("Tol0.1")){
    fTolerance = 0.1;
  }
  else if(pwd.Contains("Tol0.5")){
    fTolerance = 0.5;
  }
  else if(pwd.Contains("Tol1.0")){
    fTolerance = 1.0;
  }
  else if(pwd.Contains("Tol2.0")){
    fTolerance = 2.0;
  }
  else{
    printf("dEdxGausN::Ini Tol (fTolerance) not specified %s\n", pwd.Data()); exit(1);
  }

  //================================= fRegOpt
  fRegOpt=0;
  if(pwd.Contains("RegOpt1")){
    fRegOpt=1;
  }
  else if(pwd.Contains("RegOpt2")){
    fRegOpt=2;
  }
  else if(pwd.Contains("RegOpt3")){
    fRegOpt=3;
  }
  else if(pwd.Contains("RegOpt4")){
    fRegOpt=4;
  }

  //==================================== fkFixMean
  fkFixMean = kFALSE;
  if(pwd.Contains("FixMean")){
    fkFixMean = kTRUE;
  }

  //==================================== fElectronHighPTrigger fElectronHighPFrac
  fElectronHighPTrigger = -999;
  fElectronHighPFrac = -999;
  if(pwd.Contains("FixHighPElectron")){
    if(pwd.Contains("FixHighPElectronMode0")){
      fElectronHighPTrigger = 10;
      fElectronHighPFrac = 5e-3;
    }
    else if(pwd.Contains("FixHighPElectronMode1")){
      fElectronHighPTrigger = 10;
      fElectronHighPFrac = 1.5e-2;
    }
    else if(pwd.Contains("FixHighPElectronMode2")){
      fElectronHighPTrigger = 20;
      fElectronHighPFrac = 1.0e-2;
    }
    else{
      printf("MN2D::Ini FixHighPElectron unknown mode! %s\n", pwd.Data()); exit(1);
    }
  }

  //===================================== fRegWeight
  fRegWeight = tpchist.AnchorN();
    /*
1e9;
  if(pwd.Contains("RegWeight2E3")){
    fRegWeight = 2E3;
  }
  else if(pwd.Contains("RegWeight1E4")){
    fRegWeight = 1E4;
  }
  else{
    printf("MN2D::Ini RegWeight default %e\n", fRegWeight); 
  }
    */

  //------------------------------------------------------------ Initialize Cp initial value ----------
  fCpIni = -999;
  if(pwd.Contains("CpIniN1")){
    fCpIni = -0.01;
  }
  else if(pwd.Contains("CpIniN2")){
    fCpIni = -0.02;
  }
  else if(pwd.Contains("CpIniN3")){
    fCpIni = -0.03;
  }
  else if(pwd.Contains("CpIniN4")){
    fCpIni = -0.04;
  }
  else if(pwd.Contains("CpIniN5")){
    fCpIni = -0.05;
  }
  else if(pwd.Contains("CpIni0")){
    fCpIni = 0;
  }
  else if(pwd.Contains("CpIniP1")){
    fCpIni = 0.01;
  }
  else if(pwd.Contains("CpIniP2")){
    fCpIni = 0.02;
  }
  else{
    printf("MN2D::Ini unknown CpIni %s\n", pwd.Data()); exit(1);
  }


  //--------------------------------------------------------- physical initialization

  fHoriginal = horiginal;

  TH1D * hdummyyield = 0x0;
  fHrest = tpchist.GetFitData(fHFitData, hdummyyield, fHprefilter, fHFitYield, 1);
  delete hdummyyield;

  fHrestyield = BaseUtils::GetHyield(fHrest, "hrestyield");

  const Int_t nxrb = tpchist.Xrebin();

  //true yield instead of raw yield, the latter is of special binning and therefore with wrong normalization in low mom part
  fHTrueYield = horiginal->ProjectionX("trueyield",0, horiginal->GetNbinsY()+1);
  if(nxrb>1){
    fHTrueYield->RebinX(nxrb);
  }
  //---

  fNbin = nb;

  if(fNbin<=fNPolint){
    fNPolint = ((fNbin-1)/2)*2; //maximum Npolint
    //printf("MN2D::GetFCNChisquare fNbin<=fNPolint!! %d %d\n", fNbin, fNPolint); exit(1);
  }

  fBinX = new Int_t[fNbin];
  memcpy(fBinX, binx, sizeof(Int_t)*fNbin); 

  fXCenter = new Double_t[fNbin];
  fFitYieldArray = new Double_t[fNbin];
  fTrueYieldArray = new Double_t[fNbin];
  for(Int_t ibin=0;ibin<fNbin; ibin++){
    const Int_t ix = fBinX[ibin];

    fXCenter[ibin] = fHFitData->GetXaxis()->GetBinCenter(ix);
    fFitYieldArray[ibin] = fHFitYield->GetBinContent(ix);

    const Int_t iTrue = fHTrueYield->GetXaxis()->FindBin(fXCenter[ibin]);
    if( TMath::Abs(fHTrueYield->GetXaxis()->GetBinCenter(iTrue)-fXCenter[ibin])>EPSILON ){
      printf("MN2D::Ini wrong trueyield x ibin %d iTrue %d xx %e %e\n", ibin, iTrue, fXCenter[ibin], fHTrueYield->GetXaxis()->GetBinCenter(iTrue)); exit(1);
    }
    fTrueYieldArray[ibin] = fHTrueYield->GetBinContent(iTrue);
    printf("MN2D::Ini TrueYield ibin %d iTrue %d xx %e count %e\n", ibin, iTrue, fXCenter[ibin], fTrueYieldArray[ibin]);
  }

  //prepare preID
  fkPreID = kFALSE;
  if(hpreID){
    fkPreID = kTRUE;
    fPrePar = new Double_t[fNbin*dEdxGausN::Ntype()];
    fPreAll = new Double_t[fNbin];

    if(nxrb>1){
      hpreID->RebinX(nxrb);
    }
    
    for(Int_t ibin=0; ibin<fNbin; ibin++){
      const Int_t iPre = hpreID->GetXaxis()->FindBin(fXCenter[ibin]);
      if(TMath::Abs(hpreID->GetXaxis()->GetBinCenter(iPre)-fXCenter[ibin])>EPSILON){
        printf("MN2D::Ini hpreID X bin center not correct!! %d %d %e %e\n", ibin, iPre, hpreID->GetXaxis()->GetBinCenter(iPre), fXCenter[ibin]); exit(1);
      }

      fPreAll[ibin] = 0;

      for(Int_t itype=0; itype<dEdxGausN::Ntype(); itype++){
        const Int_t tmpy = hpreID->GetYaxis()->FindBin(itype+EPSILON);
        const Double_t ntmp = hpreID->GetBinContent(iPre, tmpy);
        fPrePar[ibin*dEdxGausN::Ntype()+itype] = ntmp;
        fPreAll[ibin] += ntmp;
        printf("MN2D::Ini preID ibin %d itype %d par %e all %e\n", ibin, itype, fPrePar[ibin*dEdxGausN::Ntype()+itype], fPreAll[ibin]);
      }
    }
  }
  //------------------------

  fYbw = fHFitData->GetYaxis()->GetBinWidth(1);

  fNYbins = fHFitData->GetNbinsY();
  fYCenter = new Double_t[fNYbins];
  for(Int_t iy = 1; iy<= fNYbins; iy++){
     fYCenter[iy] = fHFitData->GetYaxis()->GetBinCenter(iy);
  }

  fIfixed = new Bool_t[Npar()];
  SetIfixed();
 
  fLows = new Double_t[Npar()];
  fHighs= new Double_t[Npar()];
  SetLowHigh();

  fPars = new Double_t[Npar()];
  memcpy(fPars,              bbpar,  sizeof(Double_t)*dEdxGausN::NparBB());

  IniResPar(&(fPars[dEdxGausN::NparBB()]), realRes);

  memcpy(&(fPars[NparBuffer()]), binpar, sizeof(Double_t)*fNbin*NparBin());

  //must go after initialization of fPars
  FixHighPElectron();

  //------------------------------------------------------------ Initialize crossing region ---------------
  fMinXing = -999;
  fMaxXing = -999;
  SetCrossing();

  //---------------------------------------------------------------- Dummy Initialization

  fLowErrs = new Double_t[Npar()];
  fHighErrs = new Double_t[Npar()];

  for(Int_t ii=0; ii<Npar(); ii++){
    fLowErrs[ii]=0;
    fHighErrs[ii]=0;
  }

  //only correlation of A's with pion is saved!!
  fAs = new Double_t[fNbin*dEdxGausN::Ntype()];
  fErrAs = new Double_t[fNbin*dEdxGausN::Ntype()];
  fRatio = new Double_t[fNbin*dEdxGausN::Ntype()];
  fErrRatio = new Double_t[fNbin*dEdxGausN::Ntype()];
  fACov = new Double_t[fNbin*dEdxGausN::Ntype()];

  for(Int_t ii=0; ii<fNbin*dEdxGausN::Ntype(); ii++){
    fAs[ii]=0;
    fErrAs[ii]=0;
    fRatio[ii]=0;
    fErrRatio[ii]=0;
    fACov[ii]=0;
  }

  fLog = fopen("logmn2d.txt","w");
  fStopwatch = new TStopwatch;

  Print("MN2D::Ini stdout pre GetFCNChisquare",1,0,0);

  fBeta = 0.0;
  GetFCNChisquare(fPars);

  Print("MN2D::Ini fLog");
}

Double_t MN2D::SetCrossing()
{
  //the initial dEdx value should be accurate enough to determine the crossing region, which is only used by modification of npolint
 
  const Double_t xe0 = dEdxGausN::FindCrossingLog10P(fPars, dEdxGausN::BBFunc(), dEdxGausN::Mass(dEdxGausN::kELECTRON),   dEdxGausN::Mass(dEdxGausN::kKAON));
  const Double_t xe1 = dEdxGausN::FindCrossingLog10P(fPars, dEdxGausN::BBFunc(), dEdxGausN::Mass(dEdxGausN::kELECTRON),   dEdxGausN::Mass(dEdxGausN::kPROTON));
 
  const Double_t x0 = dEdxGausN::FindCrossingLog10P(fPars, dEdxGausN::BBFunc(), dEdxGausN::Mass(dEdxGausN::kPION),   dEdxGausN::Mass(dEdxGausN::kKAON));
  const Double_t x1 = dEdxGausN::FindCrossingLog10P(fPars, dEdxGausN::BBFunc(), dEdxGausN::Mass(dEdxGausN::kPROTON), dEdxGausN::Mass(dEdxGausN::kPION));
  const Double_t x2 = dEdxGausN::FindCrossingLog10P(fPars, dEdxGausN::BBFunc(), dEdxGausN::Mass(dEdxGausN::kKAON),   dEdxGausN::Mass(dEdxGausN::kPROTON));

  printf("MN2D::SetCrossing log10P %e %e, %e %e %e, mom %e %e, %e %e %e\n", xe0, xe1, x0, x1, x2, TMath::Power(10,xe0), TMath::Power(10,xe1), TMath::Power(10, x0), TMath::Power(10, x1), TMath::Power(10, x2));

  if( xe0 > xe1 || xe0 > x0 || xe0 > x1 || xe0 > x2){
    printf("MN2D::SetCrossing unexpected ordering xe0 > xe1 || xe0 > x0 || xe0 > x1 || xe0 > x2 %e %e %e %e %e\n",xe0, xe1, x0, x1, x2);
    exit(1);
  }
  if( x2 < x1 || x2 < x0 || x2 < xe0 || x2 < xe1){
    printf("MN2D::SetCrossing unexpected ordering  x2 < x1 || x2 < x0 || x2 < xe0 || x2 < xe1 %e %e %e %e %e\n", xe0, xe1, x0, x1, x2);
    exit(1);
  }
  fMinXing = xe0;
  fMaxXing = x2;
}
 
void MN2D::Clear()
{
  delete[] fBinX;
  delete[] fXCenter;
  delete[] fFitYieldArray;
  delete[] fTrueYieldArray;
  delete[] fPrePar;
  delete[] fPreAll;
  delete[] fYCenter;
  delete[] fPars;
  delete[] fAs;
  delete[] fErrAs;
  delete[] fRatio;
  delete[] fErrRatio;
  delete[] fACov;
  delete[] fLowErrs;
  delete[] fHighErrs;
  delete[] fIfixed;
  delete[] fLows;
  delete[] fHighs;

  fclose(fLog);
  //FILE can not be deleted, otherwise crashes!
  //delete fLog;

  delete fStopwatch;
  delete fMinuit;
}

//==============================================================================
//                                    functions
//==============================================================================
Double_t MN2D::IDdEdxM1Q0(const Double_t *par, const Double_t mean)
{
  return TMath::Abs(par[0]) * ( 1 + par[1]*mean + par[2]*mean*mean );
}

Double_t MN2D::IDdEdxM2Q0(const Double_t *par, const Double_t mean)
{
  //rewritting the power term in denominator and add eps to protect!
  return TMath::Abs(par[0]) * ( 1/( BaseUtils::ProtectedPow(mean,1-par[1])+EPSILON )+ par[2]/mean );
}

Double_t MN2D::IDdEdxM3Q0(const Double_t *par, const Double_t mean)
{
  return TMath::Abs(par[0]) * ( 1 + par[1]/TMath::Sqrt(mean) + par[2]/mean);
}

//==============================================================================
//                                    parameters
//==============================================================================
Int_t MN2D::NparBuffer()
{
  return dEdxGausN::NparBB()+NparRes();
}

Int_t MN2D::NparRes()
{
  //fkParaResolution = 1000, 2000, 3000
  return 3;
}

void MN2D::SetIfixed()
{
  for(Int_t ipar=0; ipar<Npar(); ipar++){
    fIfixed[ipar] = kFALSE;
  }
  fNfix = 0;

  if( dEdxGausN::Ntype()==3 && dEdxGausN::BBOpt()==dEdxGausN::kmBB){
    fIfixed[dEdxGausN::NparBB()-1] =  kTRUE;
    fNfix++;
  }

  if(fkFixMean){
    fNfix += dEdxGausN::SetIfixAllBB(fIfixed, kTRUE);
  }

  for(Int_t ibin=0; ibin<fNbin; ibin++){
    Bool_t *mgpfix=&(fIfixed[NparBuffer()+ibin*NparBin()]);

    for(Int_t itype=0; itype<dEdxGausN::Ntype(); itype++){
      mgpfix[itype*dEdxGausN::NparType()+1] = kTRUE; fNfix++;
      mgpfix[itype*dEdxGausN::NparType()+2] = kTRUE; fNfix++;
    }
  }
}

void MN2D::FixOverflow()
{
  //effect: low momentum proton and kaon should have fraction = 0
  
  Double_t *localpar = fPars;

  for(Int_t ibin=0;ibin<fNbin; ibin++){
    const Int_t ix = fBinX[ibin];
    const Double_t xx = fXCenter[ibin]; 
    
    Int_t order[100];
    const Double_t * binpar = GetBinPar(ibin, xx, localpar, order);

    //resotre to natural gaus par, order not changed!!
    Double_t gauspar[NparBin()];
    dEdxGausN::ToGausPar(binpar, gauspar);

    Double_t *mgpar = &(fPars[   NparBuffer()+ibin*NparBin()]);
    Bool_t *mgpfix  = &(fIfixed[NparBuffer()+ibin*NparBin()]);
    
    for(Int_t itype=0; itype<dEdxGausN::Ntype(); itype++){
      //has to use GausPar
      if(gauspar[itype*dEdxGausN::NparType()+1] > (TPCHist::Y1()-EPSILON)){
        mgpar[ itype*dEdxGausN::NparType()] = 0.0;
        if(!mgpfix[itype*dEdxGausN::NparType()]){
          mgpfix[itype*dEdxGausN::NparType()] = kTRUE; fNfix++;
          printf("MN2D::FixOverflow() fixing ibin %d ix %d xx %f itype %d mean %f fNfix %d\n", ibin, ix, xx, itype, gauspar[itype*dEdxGausN::NparType()+1], fNfix);
        }
      }
      else{
        if(mgpfix[itype*dEdxGausN::NparType()]){
          mgpfix[itype*dEdxGausN::NparType()] = kFALSE; fNfix--;
          printf("MN2D::FixOverflow() de-fixing ibin %d ix %d xx %f itype %d mean %f fNfix %d\n", ibin, ix, xx, itype, gauspar[itype*dEdxGausN::NparType()+1], fNfix);
        }
      }
    }
  }
}

Int_t MN2D::SetIfixAllRes(Bool_t pfix[], const Bool_t kfix)
{
  const Int_t nfix = NparRes();
  Int_t nflip = 0;
  for(Int_t ii=0; ii<nfix; ii++){
    const Int_t dd = dEdxGausN::NparBB()+ii;
    if(pfix[dd] != kfix){
      pfix[dd] = kfix;
      nflip++;
    }
  }

  if(kfix)
    return nflip;
  else
    return nflip*(-1);
}

void MN2D::FixHighPElectron()
{
  if(dEdxGausN::Ntype()==3)
    return;

  if(fElectronHighPTrigger<0)
    return;

  printf("\n");
  const Double_t triggerX = TMath::Log10(fElectronHighPTrigger);
  for(Int_t ibin = 0; ibin<fNbin; ibin++){
    const Int_t ix = fBinX[ibin];
    const Double_t xx = fXCenter[ibin];
    if(xx<triggerX)
      continue;
    
    printf("MN2D::FixHighPElectron fixing %d %d %e %e\n", ibin, ix, xx, TMath::Power(10,xx));
    fIfixed[NparBuffer()+ibin*NparBin()+dEdxGausN::kELECTRON*dEdxGausN::NparType()+0] = kTRUE; fNfix++;
    fPars[  NparBuffer()+ibin*NparBin()+dEdxGausN::kELECTRON*dEdxGausN::NparType()+0] = fElectronHighPFrac;
  }
  printf("\n");
}

void MN2D::SetLowHigh()
{
  for(Int_t ipar=0; ipar<Npar(); ipar++){
    fLows[ipar] = 0;
    fHighs[ipar]= 0;
  }

  Double_t *bslow =&(fLows[ NparBuffer()]);
  Double_t *bshigh=&(fHighs[NparBuffer()]);
  for(Int_t ibin=0; ibin<fNbin; ibin++){
    for(Int_t itype=0; itype<dEdxGausN::Ntype(); itype++){
      //---> limits will be released in nbeta>0
      bslow[ ibin*NparBin()+itype*dEdxGausN::NparType()+0] = 0;
      bshigh[ibin*NparBin()+itype*dEdxGausN::NparType()+0] = 1;     
    }
  }
}

void MN2D::IniResPar(Double_t * respar, const Double_t realRes)
{
  Double_t * tmppar = &(respar[0]);

  const Int_t dedxbit = fkParaResolution/1000;
  //dEdx{2,3}xxx
  if(dedxbit==1){
    tmppar[0]=realRes;
    tmppar[1]=0;
    tmppar[2]=0;
  }
  else if(dedxbit==2){
    tmppar[0]=realRes;
    tmppar[1]=1;
    tmppar[2]=0;
  }
  else if(dedxbit==3){
    tmppar[0]=realRes;
    tmppar[1]=0;
    tmppar[2]=0;
  }
  else if(dedxbit){
    printf("MN2D::IniResPar wrong fkParaResolution %d\n", fkParaResolution); exit(1);
  }
}

Double_t MN2D::GetResolutionParametrization(const Double_t mean, const Double_t * rawParRes)
{
  //remember to updata NparRes() and IniResPar()

  const Double_t *par = &(rawParRes[0]);

  Double_t res = -999;

  switch(fkParaResolution){
  case 1000:
    res = IDdEdxM1Q0(par, mean);
    break;
  case 2000:
    res = IDdEdxM2Q0(par, mean);
    break;
  case 3000:
    res = IDdEdxM3Q0(par, mean);
    break;
  default:
    printf("MN2D::GetResolutionParametrization wrong fkParaResolution %d\n", fkParaResolution); exit(1);
  }

  //cut out negative branched to improve minimum searching, even though negative values protected by abs in resolutionGaus
  if(res<EPSILON)
    return EPSILON;
  else 
    return res;
}

const Double_t * MN2D::GetBinPar(const Int_t ibin, const Double_t xx, Double_t *glopar, Int_t *outorder, const Bool_t kprint)
{  
  //
  //glopar is modified 
  //1) in buffer for mean from output after SetMean no matter whenther fixmean=1 -> mean in each momentum bin has to be updated using the fit parameter in the buffer
  //2) and resolution from GetResolutionParameerization

  //to speed up
  /*
  if(kprint){
    printf("getbinpar before %e %e %e %e %e\n", glopar[0], glopar[1], glopar[2], glopar[3], glopar[4]);
  }
  */

  Double_t *binpar=&(glopar[NparBuffer()+ibin*NparBin()]);

  Int_t myorder[dEdxGausN::Ntype()];
 
  //============= from buffer BB part to ibin Mean in binpar
  //put bbmean to binpar
  //put order to myorder
  //glopar (BBpar) not changed

  //to speed up
  /*
  if(kprint){
    printf("getbinpar after %e %e %e %e %e\n", glopar[0], glopar[1], glopar[2], glopar[3], glopar[4]);
  }
  */

  dEdxGausN::SetMean(xx, glopar, binpar, myorder, kprint);

  //============= resolution from buffer just after BB
  Double_t * respar = &(glopar[dEdxGausN::NparBB()]);

  //=============
  Double_t realmean = 0;
  for(Int_t ith=0; ith<dEdxGausN::Ntype(); ith++){
    const Int_t gausanchor = ith*dEdxGausN::NparType();

    //>>>>>>>>>>>>>>>>> calc. width and set width to glopar(i.e. binpar) <<<<<<<<<<<<<<<<<<<<<<<
    //big bug fixed!!
    realmean += binpar[gausanchor+1];
    
    //input is itype, i.e. myorder[ith]=itype
    const Double_t calcRes = GetResolutionParametrization(realmean, respar);
    binpar[gausanchor+2] = calcRes;    
  }

  if(outorder){
    memcpy(outorder,myorder, sizeof(myorder));
  }

  //to speed up
  /*
  if(kprint){
    for(Int_t ith=0; ith<dEdxGausN::Ntype(); ith++){
      for(Int_t ip=0; ip<dEdxGausN::NparType(); ip++){
        const Int_t ipar = ith*dEdxGausN::NparType()+ip;
        printf("ibin %d xx %e ith %d ip %d ipar %d -- %e\n", ibin, xx, ith, ip, ipar, binpar[ipar]);
      }
    }
  }
  */
  return binpar;
}

void MN2D::SetPhysics(const TMatrixDSym * fullcov, const Bool_t khesse)
{
  //in each pbin, do
  //1. get the order, find where is pion
  //2. set As and ErrAs, should be identical to the one in ::Save
  //3. calc. ratio over pion
  //4. calc. err of ratio from fullcov
  //5. calc. sumcov

  for(Int_t ibin = 0; ibin < fNbin; ibin++){
    const Int_t ix = fBinX[ibin];
    const Double_t xx = fXCenter[ibin];

    Double_t dummy[NparBin()];//big bug fixed!!
    //itype = order[ith]
    Int_t order[dEdxGausN::Ntype()];
    dEdxGausN::SetMean(xx, fPars, dummy, order);

    //ith=invmap[itype]
    Int_t invmap[dEdxGausN::Ntype()];
    dEdxGausN::GetInvMap(order, invmap);

    Double_t sumcov=0;
    for(Int_t ith=0; ith<dEdxGausN::Ntype(); ith++){
      const Int_t isave = ibin*dEdxGausN::Ntype()+order[ith];
      const Int_t ipar  = NparBuffer() + ibin*NparBin() + ith*dEdxGausN::NparType();
      //printf("longtest %d %d -- %d %d -- %d\n", ith, order[ith], ipar, ibin*dEdxGausN::Ntype()+order[ith], ibin);

      fAs[isave] = TMath::Abs(fPars[ipar]);
     
      if(khesse){

        //use binomial error for fraction, NOT using TMinuit error to avoid inflation by NPOLINT>2, but the TMinuit covariance is still needed for ratio error !!
        Double_t tmpeff =  fAs[isave];
        if( fAs[isave] <EPSILON ){
          tmpeff = EPSILON;
        }
        if( fAs[isave] > 1-EPSILON ){
          tmpeff = 1-EPSILON;
        }
        
        fErrAs[isave] =  TMath::Sqrt(tmpeff*(1-tmpeff)/fTrueYieldArray[ibin]);

        /*
      const Double_t sigma2 = (*fullcov)[ipar][ipar];
      if(sigma2<EPSILON){
        if(TMath::Abs(sigma2)<EPSILON){
          fErrAs[isave] = 0;
        }
        else{
          printf("MN2D::SetPhysics sigma2<EPSILON %e %d\n", sigma2, ipar); exit(1);
        }
      }
      else{
        //fErrAs[isave] = TMath::Sqrt(sigma2);
      }
        */
      }

      const Int_t ipion = NparBuffer() + ibin*NparBin() + invmap[dEdxGausN::kPION]*dEdxGausN::NparType();
      fACov[isave] = (*fullcov)[ipar][ipion];

      //sum of cov elements and then sqrt is error fo \sum A
      for(Int_t jth=0; jth<dEdxGausN::Ntype(); jth++){
        const Int_t jpar   = NparBuffer() + ibin*NparBin() + jth*dEdxGausN::NparType();
        sumcov+=(*fullcov)[ipar][jpar];
      }
    }

    for(Int_t itype=0; itype<dEdxGausN::Ntype(); itype++){
      const Int_t ia = ibin*dEdxGausN::Ntype()+itype;
      const Int_t ib = ibin*dEdxGausN::Ntype()+dEdxGausN::kPION;
      const Double_t aa = fAs[ia];
      const Double_t bb = fAs[ib];
      if(bb<EPSILON){
        fRatio[ia]=fErrRatio[ia]=-999;
      }
      else{
        const Double_t ff = aa/bb;
        fRatio[ia]=ff;

        if(aa>EPSILON){
          const Double_t v0 = ff/aa;
          const Double_t v1 = -ff/bb;
          fErrRatio[ia]=BaseUtils::GetCorrelatedError(v0,v1,TMath::Power(fErrAs[ia],2), fACov[ia], TMath::Power(fErrAs[ib],2));
        }
        else{
          fErrRatio[ia]=0;
        }
      }
    }
  }
}

//==============================================================================
//                                    output
//==============================================================================
Char_t* MN2D::Status(Char_t * info, const Bool_t kend)
{
  Char_t * status = Form("fkErr = %d, fStatChi2 = %e, NDF() = %d, ReducedChi2() = %e, fBeta =  %e, fPenalty = %e, b*p = %e fcn = %e", fkErr, fStatChi2, NDF(), ReducedChi2(), fBeta, fPenalty, fBeta*fPenalty, myFCN());
  printf("%-110s %s%c", info, status, kend?'\n':' ');
  return status;
}

void MN2D::Print(const TString tag, const Bool_t kprintBinPar, const Bool_t kprintBinErr, const Bool_t ktofile)
{
  FILE *pout = stdout;
  if(ktofile)
    pout=fLog;

  fprintf(pout, "\n================================================ %s:\n", tag.Data());

  fprintf(pout, "fkParaResolution = %d;\n", fkParaResolution);
  fprintf(pout, "fNPolint = %d;\n", fNPolint);
  fprintf(pout, "fTolerance = %e;\n", fTolerance);

  fprintf(pout, "fRegOpt = %d;\n", fRegOpt);
  fprintf(pout, "fkFixMean = %d;\n", fkFixMean);

  fprintf(pout, "fElectronHighPTrigger = %e\n", fElectronHighPTrigger);
  fprintf(pout, "fElectronHighPFrac = %e\n", fElectronHighPFrac);

  fprintf(pout, "fRegWeight = %e\n", fRegWeight);
  fprintf(pout, "fCpIni = %e\n", fCpIni);
  fprintf(pout, "fkUseTrueYield = %d\n", fkUseTrueYield);
  
  fprintf(pout, "NparBuffer(): %d;\n", NparBuffer());
  fprintf(pout, "dEdxGausN::NparBB(): %d;\n", dEdxGausN::NparBB());
  fprintf(pout, "NparRes(): %d\n", NparRes());

  fprintf(pout, "NparBin(): %d;\n", NparBin());
  fprintf(pout, "Npar(): %d;\n", Npar());
  fprintf(pout, "Nfree(): %d\n", Nfree());
  fprintf(pout, "fNfix = %d\n", fNfix);

  fprintf(pout, "fkPreID = %d\n", fkPreID);
  fprintf(pout, "RegularizationWeight(1e6) %e\n", RegularizationWeight(1e6));
  fprintf(pout, "RegMaxdEdx() %e\n", RegMaxdEdx());

  fprintf(pout, "fNdata = %d\n", fNdata);
  fprintf(pout, "Status: %s\n", Status(Form("-> MN2D::Print %s", tag.Data())));

  fprintf(pout, "fNbin = %d;\n", fNbin);

  fprintf(pout, "fMinXing = %e;\n", fMinXing);
  fprintf(pout, "fMaxXing = %e;\n", fMaxXing);

  for(Int_t ii=0; ii<NparBuffer(); ii++){
    fprintf(pout, "fPars[%d] = %e;\n", ii, fPars[ii]);
  }
  
  fprintf(pout, "\n");

  if(!kprintBinPar)
    return;

  for(Int_t ii=0; ii<fNbin; ii++){
    fprintf(pout, "fBinX[%3d] = %d; ", ii, fBinX[ii]);
    for(Int_t ip = 0; ip<NparBin(); ip++){
      const Int_t id = NparBuffer()+ii*NparBin()+ip;
      fprintf(pout, "fPars[%3d] = %10e; ", id, fPars[id]);
    }
    fprintf(pout, "\n");
  }

  fprintf(pout, "\n");

  if(!kprintBinErr)
    return;

  for(Int_t ii=0; ii<fNbin; ii++){
    fprintf(pout, "fBinX[%3d] = %d; ", ii, fBinX[ii]);
    for(Int_t ip = 0; ip<NparBin(); ip++){
      const Int_t id = NparBuffer()+ii*NparBin()+ip;
      fprintf(pout, "fLowErrs[%3d] = %10e; ", id, fLowErrs[id]);
    }
    fprintf(pout, "\n");
  }

  fprintf(pout, "\n");

  for(Int_t ii=0; ii<fNbin; ii++){
    fprintf(pout, "fBinX[%3d] = %d; ", ii, fBinX[ii]);
    for(Int_t ip = 0; ip<NparBin(); ip++){
      const Int_t id = NparBuffer()+ii*NparBin()+ip;
      fprintf(pout, "fHighErrs[%3d] = %10e; ", id, fHighErrs[id]);
    }
    fprintf(pout, "\n");
  }

  fprintf(pout, "\n");

  PrintPhysics(pout, "fAs", fAs);
  PrintPhysics(pout, "fErrAs", fErrAs);
  PrintPhysics(pout, "fACov", fACov);
  PrintPhysics(pout, "fRatio", fRatio);
  PrintPhysics(pout, "fErrRatio", fErrRatio);  

  fprintf(pout, "\n");
}

void MN2D::PrintPhysics(FILE *pout, const TString ss, const Double_t * var)
{
  for(Int_t ibin=0; ibin<fNbin; ibin++){
    fprintf(pout, "bin=%4d; ", ibin);
    for(Int_t itype=0; itype<dEdxGausN::Ntype(); itype++){
      const Int_t id = ibin*dEdxGausN::Ntype()+itype;
      fprintf(pout, "%s[%3d] = %10e; ", ss.Data(), id, var[id]);
    }
    fprintf(pout, "\n");
  }
  fprintf(pout, "\n");
}

void MN2D::Save(const TString tag)
{
  TTreeSRedirector *fStream = new TTreeSRedirector(Form("%s_save.root", tag.Data()));

  for(Int_t ip=0; ip<NparBuffer(); ip++){
    (*fStream)<<"tree"<<
      Form("bb%d=",ip)<<fPars[ip]<<
      Form("leb%d=",ip)<<fLowErrs[ip]<<
      Form("heb%d=",ip)<<fHighErrs[ip];
  }
  Double_t chidof = ReducedChi2();
  Int_t npar = Npar();
  Int_t ndf = NDF();
  (*fStream)<<"tree"<<
    "npar="<<npar<<
    "chi2="<<fStatChi2<<
    "ndof="<<ndf<<
    "chidof="<<chidof<<
    "penalty="<<fPenalty<<
    "nbin="<<fNbin<<
    "\n";

  //localpar, i.e. fPars will be updated in GetBinPar (bufferResPart, mean and resolution are set).
  Double_t *localpar = fPars;
  Double_t *locallowerr  = fLowErrs;
  Double_t *localhigherr = fHighErrs;

  for(Int_t ibin=0; ibin<fNbin; ibin++){
    Int_t ix = fBinX[ibin];
    Double_t xx = fXCenter[ibin];

    //need to put order here, otherwise the output has wrong order
    //only record order, no reordering
    Int_t order[dEdxGausN::Ntype()];
    GetBinPar(ibin, xx, localpar, order);

    //totalyield, outpar, outlowerr, outhigherr have to be defined here due to streamer
    Double_t yields[dEdxGausN::Ntype()];
    Double_t outpar[npar];
    Double_t outlowerr[npar];
    Double_t outhigherr[npar];
    Double_t totalyield = fTrueYieldArray[ibin];
    const Int_t ipar = NparBuffer()+ibin*NparBin();
    dEdxGausN::SavePar(fStream, "mn2d", &totalyield, yields, &(localpar[ipar]), &(locallowerr[ipar]), &(localhigherr[ipar]), order, outpar, outlowerr, outhigherr);

    SavePhysics(fStream,"mn2d", ibin);

    (*fStream)<<"mn2d"<<
      "ix="<<ix<<
      "xx="<<xx<<   
      "\n";
  }
  
  fStream->GetFile()->cd();
  fHFitData->Write();
  delete fStream;

  Print(Form("MN2D::Save %s", tag.Data()));
}

void MN2D::SavePhysics(TTreeSRedirector * stream, const TString tn, const Int_t ibin)
{
  for(Int_t itype=0; itype<dEdxGausN::Ntype(); itype++){
    const Int_t id = ibin*dEdxGausN::Ntype()+itype;
    (*stream)<<tn<<
      Form("aa%d=",   itype)<<fAs[id]<<
      Form("ea%d=",   itype)<<fErrAs[id]<<
      Form("cova%d=", itype)<<fACov[id]<<
      Form("rr%d=",   itype)<<fRatio[id]<<
      Form("er%d=",   itype)<<fErrRatio[id];
  }
}

void MN2D::Draw2d(const TString tag)
{
  //_______________________________________ define output _____________________
  TList *lkept = new TList;
  lkept->SetOwner(kFALSE);

  lkept->Add(fHoriginal);
  lkept->Add(fHprefilter);
  lkept->Add(fHFitData);
  lkept->Add(fHFitYield);
  lkept->Add(fHTrueYield);
  lkept->Add(fHrest);
  lkept->Add(fHrestyield);

  TList *ll = new TList;
  ll->SetOwner(kTRUE);

  TH2D *hfit = (TH2D*)fHFitData->Clone(tag+"_fit");
  hfit->Scale(0);
  ll->Add(hfit);

  const Double_t * fitxarr= BaseUtils::GetAxisArray(hfit->GetXaxis());
  if(!fitxarr){
    printf("MN2D::Draw2d null fitxarr %d!!\n", hfit->GetNbinsX()); exit(1);
  }

  TH1D *hmean[dEdxGausN::Ntype()];
  for(Int_t itype = 0; itype<dEdxGausN::Ntype(); itype++){
    hmean[itype] = new TH1D(tag+Form("_mean%d", itype),"",hfit->GetNbinsX(), hfit->GetXaxis()->GetXmin(), hfit->GetXaxis()->GetXmax());
    //big bug fixed!
    hmean[itype]->GetXaxis()->Set(hfit->GetNbinsX(), fitxarr);
 
    //hmean[itype]->GetXaxis()->SetRange(fBinX[0], fBinX[fNbin-1]);
    ll->Add(hmean[itype]);
  }

  TH1D * hchi2 = new TH1D(tag+"_chi2","",hfit->GetNbinsX(), hfit->GetXaxis()->GetXmin(), hfit->GetXaxis()->GetXmax());
  //big bug fixed!
  hchi2->GetXaxis()->Set(hfit->GetNbinsX(), fitxarr);

  hchi2->GetXaxis()->SetRange(fBinX[0], fBinX[fNbin-1]);
  hchi2->SetMarkerStyle(20);
  ll->Add(hchi2);

  //______________________________________________________________________
  //______________________________________________________________________

  //localpar is updated in GetBinPar, (bufferResPart, mean and resolution is set).
  Double_t *localpar = fPars;

  for(Int_t ibin=0;ibin<fNbin; ibin++){
    const Int_t ix = fBinX[ibin];
    const Double_t xx = fXCenter[ibin]; 

    Int_t order[100];
    const Double_t * binpar = GetBinPar(ibin, xx, localpar, order);

    //resotre to natural gaus par, order not changed!!
    Double_t gauspar[NparBin()];
    dEdxGausN::ToGausPar(binpar, gauspar);
    Int_t invmap[100];
    dEdxGausN::GetInvMap(order, invmap);
    
    Double_t chi2bini = 0;
    Double_t ndfbini = 0;
    for(Int_t iy = 1; iy<= hfit->GetNbinsY(); iy++){
      const Double_t yy = hfit->GetYaxis()->GetBinCenter(iy);
      const Double_t cc = dEdxGausN::Func(&yy, binpar);

      //set hfit
      hfit->SetBinContent(ix, iy, cc); hfit->SetBinError(ix,iy,0);

      //calc chi2
      const Double_t rawc = fHFitData->GetBinContent(ix, iy);
      const Double_t rawerr = fHFitData->GetBinError(ix, iy);
      if(rawerr>EPSILON){
        chi2bini+=TMath::Power((rawc-cc)/rawerr,2);
        ndfbini++;
      }

      //set mean
      for(Int_t itype = 0; itype<dEdxGausN::Ntype(); itype++){
        hmean[itype]->SetBinContent(ix, gauspar[invmap[itype]*dEdxGausN::NparType()+1]); hmean[itype]->SetBinError(ix, 0);
      }
    }

    //set chi2
    //just need a simple chi2/ndof
    //ndfbini -= Nfree()/fNbin  ;
    hchi2->SetBinContent(ix, ndfbini>EPSILON? chi2bini/ndfbini : 0); hchi2->SetBinError(ix, 0);

    //get slice
    TH1D *tmph = fHFitData->ProjectionY(Form("%sx%d", fHFitData->GetName(), ix),ix,ix,"oe");
    tmph->SetTitle(Form("x %e",xx));
    tmph->SetLineColor(kBlack);
    tmph->SetMarkerColor(kBlack);
    ll->Add(tmph);

    //draw slice
    TH1D *hslicefit = BaseUtils::GetHfit(Form("%sfit", tmph->GetName()), dEdxGausN::Func, binpar, tmph->GetXaxis()->GetXmin(), tmph->GetXaxis()->GetXmax());
    hslicefit->SetLineColor(kBlack);
    //hslicefit->SetDrawOption("same hist C");
    ll->Add(hslicefit);

    Int_t col[]={kBlue, kRed, kGreen+3, kMagenta, kYellow};
    for(Int_t itype=0; itype<dEdxGausN::Ntype(); itype++){
      dEdxGausN::DrawTF1(itype, tmph->GetName(), &(gauspar[invmap[itype]*dEdxGausN::NparType()]), col[itype], ll, fLog); 
    }
  }

  //remove empty bins
  TGraph * grmean[dEdxGausN::Ntype()];
  for(Int_t itype = 0; itype<dEdxGausN::Ntype(); itype++){
    grmean[itype] = new TGraph;
    grmean[itype]->SetName(Form("%sgr",hmean[itype]->GetName()));
    ll->Add(grmean[itype]);
    Int_t ig = 0;
    for(Int_t ii=1; ii<=hmean[itype]->GetNbinsX(); ii++){
      const Double_t yy = hmean[itype]->GetBinContent(ii);
      if(yy<EPSILON){
        continue;
      }

      const Double_t xx = hmean[itype]->GetBinCenter(ii);

      grmean[itype]->SetPoint(ig, xx, yy);
      ig++;
    }
  }

  TGraph *grchi2 = new TGraph;
  grchi2->SetName(Form("%sgr",hchi2->GetName()));
  ll->Add(grchi2);
  Int_t ig = 0;
  for(Int_t ii=1; ii<=hchi2->GetNbinsX(); ii++){
    const Double_t yy = hchi2->GetBinContent(ii);
    if(yy<EPSILON){
      continue;
    }
    
    const Double_t xx = hchi2->GetBinCenter(ii);
    
    grchi2->SetPoint(ig, xx, yy);
    printf("MN2D::Draw2d %s chi2 ig %d xx %e chi2/ndf %e\n", tag.Data(), ig, xx, yy);
    ig++;
  }

  TFile *fout= new TFile(Form("%s_draw2d.root", tag.Data()),"recreate");
  lkept->Write();
  ll->Write();
  fout->Save();
  fout->Close();
  delete fout;
  delete ll;
  delete lkept;

  Print(Form("MN2D::Draw2d %s", tag.Data()));
}

//==============================================================================
//                                    fit
//==============================================================================
Double_t MN2D::LinearFit(const Int_t np, const Double_t xarr[], const Double_t yarr[], const Double_t eys[], const Double_t xtarget, const Int_t order)
{
  //chi2 polinomial fit of degree, i.e. up to x^(deg-1), see also Recipes p673.
  const Int_t deg = order+1;
  TMatrixD leftmatrix(deg, deg);
  TMatrixD rightvector(deg,1);

  //loop over data point
  for(Int_t ip=0; ip<np; ip++){
    const Double_t xx = xarr[ip]-xtarget;
    const Double_t yy = yarr[ip];
    //to avoid ey<0 or ey=0
    const Double_t ey = TMath::Max(EPSILON, eys[ip]);

    //starting from [0][0] term
    Double_t lastrecord = 1/ey/ey;    
    //loop over each matrix entry once with increasing power of xx
    for(Int_t sumij=0; sumij<=(deg-1)*2; sumij++){
      for(Int_t ii=TMath::Max( 0, sumij-(deg-1) ); ii<=TMath::Min(deg-1, sumij); ii++){
        const Int_t jj = sumij - ii;
        //printf("%d: %d %d\n", sumij, ii, jj);
        
        leftmatrix[ii][jj] += lastrecord;
      }
      //printf("sum %d\n\n",sumij);
      lastrecord *= xx;
    }

    //loop over vector entry
    lastrecord = yy/ey/ey;
    for(Int_t ii=0; ii<deg; ii++){
      rightvector[ii][0] +=  lastrecord;
      lastrecord *= xx;
    }

    //checked by hand, correct!
    //leftmatrix.Print();
    //printf("===\n");
    //rightvector.Print();
  }

  /*
TMatrixD a(2,2)
a[0][0]=cos(3)
a[0][1]=sin(3)
a[1][0]=-sin(3)
a[1][1]=cos(3)
a.Print()
TDecompLU lu(a)
TMatrixD c=lu.Invert()
c.Print()
a.Print()
  */

  //leftmatrix
  TDecompSVD mdec(leftmatrix);
  if(!mdec.Decompose()){
    printf("MN2D::LinearFit matrix singular?\n");
    for(Int_t ii=0; ii<np; ii++){
      printf("%d: %e %e %e - %e\n", ii, xarr[ii], yarr[ii], eys[ii], xtarget);
    }
    leftmatrix.Print();
    exit(1);
  }
  const TMatrixD inv=mdec.Invert();

  Double_t par0 = 0;
  for(Int_t ii=0; ii<deg; ii++){
    par0 += inv[0][ii]*rightvector[ii][0];
  }

  //compared to LinearFitter::Eval, consistent
  return par0;
}

void MN2D::UpdateFractionWithPreID(Double_t &p0y, const Double_t tpcall, const Double_t prepar, const Double_t preall)
{
  p0y = (p0y*tpcall+prepar)/(tpcall+preall);
}

Double_t MN2D::RegularizationWeight(const Double_t ncount)
{
  //500 is too small but about to be OK, so a factor of 2
  return TMath::Min(fRegWeight, ncount);
}

Double_t MN2D::GetFCNChisquare(Double_t localpar[])
{
  //========================================================================
  //                    cal statchi2+beta*fluc
  //========================================================================

  //to speed up
  //static Int_t istep = 0;

  /*
  const Int_t kprint = 0;//(istep>= 112000 && istep < 113000);// 0;//TMath::Max(istep-46649,0);// 
  if(kprint) {
    for(int ii=0; ii<Npar(); ii++){
      printf("istep %d %d: %e %e\n", istep, ii, fPars[ii], localpar[ii]);
    }
  }
  */

  //statistical chi2, i.e. the standard one
  fStatChi2 = 0;

  fNdata = 0;

  //======================================================================
  //calculate statistical Chi2
  //and record fraction for particle types and position

  Double_t frac[dEdxGausN::Ntype()][fNbin];
  Bool_t kreg[dEdxGausN::Ntype()][fNbin];

  for(Int_t ibin=0;ibin<fNbin; ibin++){
    const Int_t ix = fBinX[ibin];
    const Double_t xx = fXCenter[ibin];

    Int_t outorder[dEdxGausN::Ntype()];
    //order[ith]=itype
    const Double_t *binpar = GetBinPar(ibin, xx, localpar, outorder /*, kprint*/);

    //calc chi2
    //still use the regulated stat.
    const Double_t norm = fkUseTrueYield? fTrueYieldArray[ibin] : fFitYieldArray[ibin];
    const Double_t ybw = fYbw; 

    //also use empty bins
    for(Int_t iy = 1; iy<= fNYbins; iy++){
      const Double_t yy = fYCenter[iy]; 

      //no change of binpar in ::Func
      const Double_t zfit = dEdxGausN::Func(&yy, binpar);

      const Double_t zdata = fHFitData->GetBinContent(ix, iy);

      //check normalization consistency ->
      /* //to speed up
      const Double_t err = fHFitData->GetBinError(ix, iy);
      if(err>EPSILON){
        //y = kN, e = k sqrt(N), pow(y/e,2) = N, 1/k = N/y should be norm*ybw
        const Double_t invk = TMath::Power(zdata/err,2)/zdata;
        if(fabs(invk/ybw-norm)>1){
          printf("MN2D::GetFCNChisquare normalization not consistent! norm %e ybw %e zdata %e err %e invk %e delta %e\n", norm, ybw, zdata, err, invk, fabs(invk/ybw-norm));exit(1);
        }
      }
      //<-
      */

      //zdata and zfit are both pdf, fraction = pdf*bin_width
      const Double_t prob = TMath::Poisson(zdata*ybw*norm, zfit*ybw*norm);

      //try to have Binomial error and better constraint on zfit (0<,<1)
      //but seems the constraint is too strong and no all parameter space is probed (<0, >1)
      //in the end the convergence, or reduced chi2, is very bad
      //so don't use it!
      //const Double_t prob = BaseUtils::BinomialD(zdata*norm, norm, zfit);

      fStatChi2 += -TMath::Log(prob>EPSILON ? prob : EPSILON);

      fNdata ++;
    }

    //to speed up
    /*
    if(kprint){
      printf("istep %d binpar %e %e %e %e %e %e outorder %d %d %d %d\n", istep, binpar[0],binpar[1],binpar[2],binpar[3],binpar[4],binpar[5],outorder[0],outorder[1],outorder[2],outorder[3]);
    }
    */

    //save frac used only in regularization! First patch then regulate. Idea: B. Hess.
    for(Int_t ith=0; ith<dEdxGausN::Ntype(); ith++){
      const Int_t itype = outorder[ith];
      Double_t tmpfrac =  TMath::Abs(binpar[ith*dEdxGausN::NparType()+0]);      //big bug fixed !! only abs() matters
      if(fkPreID){
        const Double_t tpcall = fTrueYieldArray[ibin];
        const Double_t prepar = fPrePar[ibin*dEdxGausN::Ntype()+itype];
        const Double_t preall = fPreAll[ibin];
        UpdateFractionWithPreID(tmpfrac, tpcall, prepar, preall);

        //checked with numbers from drawpar.C and is fine
        /*
        if(xx>0){
          static int bbkk = 0;
          printf("only for test ibin %d xx %e itype %d old %e tpcall %e TOFpar %e TOFall %e new %e\n", ibin, xx, itype, TMath::Abs(binpar[ith*dEdxGausN::NparType()+0]), tpcall, prepar, preall,  tmpfrac);
          bbkk++; 
          if(bbkk==12){
            exit(1);
          }
        }
        */
      }

      frac[itype][ibin]= tmpfrac;

      //if(kprint) {
      //  printf("%d frac %d %d: %d %e %e\n", istep, ibin, ith, outorder[ith], frac[outorder[ith]][ibin], binpar[ith*dEdxGausN::NparType()+0]);
      //}
    }        

    Double_t realmean=0;
    for(Int_t ith=0; ith<dEdxGausN::Ntype(); ith++){ 
      const Int_t itype = outorder[ith];
      const Int_t gausanchor = ith*dEdxGausN::NparType();  
      realmean += binpar[gausanchor+1];

      kreg[itype][ibin] = (realmean<RegMaxdEdx());
    }

  }

  //==========================================================================
  //-----------penaltiy term on fractions ----------
  //regularized on all particles, otherwise bad behavior of the free one
 
  fPenalty = 0;

  //fNPolint = 0 no penalty, checked: consistent with before
  //need to be calculated even if beta=0 in order to evaluate next beta
  //must run even fBeta=0 to get the intial penalty
  if(fNPolint){
    //for NPolint = 2
    //* X *
    //  * X *

    //for NPolint = 4
    //* * X * *, 4-point interpolation
    //0, 1, 2, ..., n-3, n-2, n-1
    //    * * X * *

    //----> fine to ues ibin++
    //for(Int_t ibin = fNPolint/2; ibin <= fNbin-1-fNPolint/2; ibin++){

    //better boundary condition
   
    for(Int_t ibin = 0; ibin <= fNbin-1; ibin++){
      //==========================================================
      //first prepare neighbouring points (X,Y)
      
      const Double_t targetX = fXCenter[ibin];
      Double_t Xpol[fNPolint];

      //better boundary condition; boundary condition has to be symmetric!!! otherwise the end points is heavily determined by extrapolation!

      Int_t tmpNpol = 2; //fNPolint;       //use npol for trival region to have least effect on the shape
      Int_t tmpRegOpt = 1;

      //modify the range npol = NPolint in the region of crossings -->
      //should be well within boundaries
      if(fXCenter[ibin]> fMinXing && fXCenter[ibin] < fMaxXing){
        tmpNpol = fNPolint;
        tmpRegOpt = fRegOpt;
      }
      //<---
     
      if(ibin< tmpNpol/2)
        tmpNpol = ibin*2;
      else if( (fNbin-1-ibin) < tmpNpol/2 )
        tmpNpol = (fNbin-1-ibin)*2;

      //trueNpol=0,2,4,...; 0 has to be skipped, relevant for all end points
      if(tmpNpol<2)
        continue;

     
      const Int_t trueNpol = tmpNpol;
      const Int_t trueRegOpt = tmpRegOpt;

      const Int_t binstart = ibin-trueNpol/2;
      /* //too speed up
      if(binstart<0){
        printf("MN2D::GetFCNChisquare binstart<0!! %d %d %d\n", binstart, ibin, trueNpol); exit(1);
      }
      */
      const Int_t binstop = ibin+trueNpol/2;
      /*//to speed up
      if(binstop>fNbin-1){
        printf("MN2D::GetFCNChisquare binstop>fNbin-1!! %d %d %d %d\n", binstart, ibin, trueNpol, fNbin); exit(1);
      }
      */

      /*
      //test
      if(istep==0){
        printf("testbinstartbinstop fNbin %d fNPolint %d trueNpol %d ibin %d binstart %d binstop %d\n", fNbin, fNPolint, trueNpol, ibin, binstart, binstop);
      }
      */

      tmpNpol = 0;
      for(Int_t ip = binstart; ip<=binstop; ip++){
        if(ip==ibin)
          continue;

        Xpol[tmpNpol] = fXCenter[ip];
        tmpNpol++;
      }

      /*
      for(Int_t ip=0; ip<fNPolint/2; ip++){
        ///think about ip=0
        Xpol[ip]           = fXCenter[ibin-fNPolint/2+ip];
        Xpol[fNPolint-1-ip]= fXCenter[ibin+fNPolint/2-ip];
      }
      */
      const Int_t ix = fBinX[ibin];

      //using FitYield too rigid where count is high ~ 1e4
      //always use fityield for regularization
      const Double_t reduYield =  fkUseTrueYield? fTrueYieldArray[ibin] : fFitYieldArray[ibin]; //RegularizationWeight( fTrueYieldArray[ibin] + (fkPreID? fPreAll[ibin] : 0)  );

      for(Int_t itype = 0; itype<dEdxGausN::Ntype(); itype++){
        if(!kreg[itype][ibin]){
          continue;
        }

        const Double_t targetY = frac[itype][ibin];
       
        Double_t Ypol[fNPolint];
        Double_t EyPol[fNPolint];

        //better boundary condition
        tmpNpol = 0;
        for(Int_t ip = binstart; ip<=binstop; ip++){
          if(ip==ibin)
            continue;

          Ypol[tmpNpol] = frac[itype][ip];
          EyPol[tmpNpol] = TMath::Sqrt(Ypol[tmpNpol]/reduYield);
          tmpNpol++;
        }

        /*
        for(Int_t ip=0; ip<fNPolint/2; ip++){
          Ypol[ip]            = frac[itype][ibin-fNPolint/2+ip];
          Ypol[fNPolint-1-ip] = frac[itype][ibin+fNPolint/2-ip];
          //if(kprint) {
          //  printf("longtest1 %d -- %d %d %d %d %d -- %e %e %e %e\n", itype, fNPolint, ip, ibin-fNPolint/2+ip, fNPolint-1-ip, ibin+fNPolint/2-ip, Ypol[ip]  , Ypol[fNPolint-1-ip], frac[itype][ibin+fNPolint/2-ip],  frac[itype][ibin-fNPolint/2+ip]);
          //}
        }

        for(Int_t ip=0; ip<fNPolint; ip++){
          EyPol[ip] = TMath::Sqrt(Ypol[ip]/reduYield);
        }
        */

        //==========================================================
        //get the interpolation from neighbours

        //it should be very fast
        //static void Polint(const Int_t n, const Double_t *xa, const Double_t *ya, const Double_t x, Double_t &y, Double_t &dy, const Int_t kprint=0);

        Double_t Yeva = -999;
        if(trueRegOpt==0){
          Double_t dummy = -999;
          BaseUtils::Polint(trueNpol, Xpol, Ypol, targetX, Yeva, dummy);
        }
        else{
          //use pol2 instead, based on Taylor expansion, suggested by M. Ivanov. Stronger regularization
          Yeva = LinearFit(trueNpol, Xpol, Ypol, EyPol, targetX, trueRegOpt);
          //to speed up
          /*
          if(kprint){
            for(Int_t ipr=0; ipr<fNPolint; ipr++){
              printf("istep %d: %d %e %e %e -- %e %e\n", istep, ipr, Xpol[ipr], Ypol[ipr], EyPol[ipr], targetX, Yeva);
            }
          }
          */
        }

        //======================================================
        //chi^2 form of penalty

        const Double_t delta = Yeva-targetY;
        //the Abs() is important, otherwise there is "nan" and not converge
        //use (targetY+Yeva)/2 instead of targetY for <frac>
        const Double_t denom   = TMath::Sqrt(TMath::Abs(targetY+Yeva)/2/reduYield); //= sqrt(AN)/N

        if(denom<EPSILON*1e-6){
          /*
          for(Int_t ibad=0; ibad<fNbin; ibad++){
            printf("bad %d %e -- %e %e %e %e\n", ibad, fXCenter[ibad], frac[0][ibad], frac[1][ibad], frac[2][ibad], frac[3][ibad]);
          }
          printf("MN2D::GetFCNChisquare denom 0 !! ib %d itype %d tX %e X %e %e %e %e itype %d Y %e %e %e %e -- tY %e Yeva %e delta %e denom %e fPenalty %e\n", ibin, itype, targetX, Xpol[0], Xpol[1], Xpol[2], Xpol[3], itype, Ypol[0], Ypol[1], Ypol[2], Ypol[3], targetY, Yeva, delta, denom, fPenalty); //exit(1);
          */
        }
        else{
          //1/2 to restore to Gaussian p.d.f, note the sigma is considered as parameter-independent
          fPenalty += TMath::Power(delta/denom, 2) / 2. ;
        }

        //if(kprint ){
        //  printf("longtest %d ib %d tX %e X %e %e %e %e itype %d Y %e %e %e %e -- tY %e Yeva %e delta %e denom %e del/den %e w2 %e fPenalty %e -- redy %e\n", istep,  ibin, targetX, Xpol[0], Xpol[1], Xpol[2], Xpol[3], itype, Ypol[0], Ypol[1], Ypol[2], Ypol[3], targetY, Yeva, delta, denom, delta/denom, TMath::Power(delta/denom,2),fPenalty, reduYield); 
        //}

        //to speed up
        /*
        if(kprint>3){
          exit(1);
        }
        */
      }
    }
  }

  //to speed up
  /*
  if(istep%1000==0){
    printf("istep %d myFCN %e pars %e %e %e %e %e -- %e %e %e\n", istep, myFCN(), localpar[0], localpar[1], localpar[2], localpar[3], localpar[4], localpar[5], localpar[6], localpar[7]);
  }
  */
  //istep++;

  return myFCN();
}

void MN2D::MinFCNChisquare(Int_t &npar, Double_t *gin, Double_t &f, Double_t *pars, Int_t iflag)
{
  f = GetFCNChisquare(pars);
}

Int_t MN2D::DoMigrad(const Int_t nbeta, const Bool_t klong)
{
  //=================================================================================
  //             minimization kernel of the coherent fit
  //             called by
  //             ::BetaFit
  //             fit in flow of "strategy" improves convergence a lot
  //             calling
  //             ::MinFCNChisquare, which directly calls 
  //             ::GetFCNChisquare
  //=================================================================================

  /*
http://wwwasdoc.web.cern.ch/wwwasdoc/WWW/minuit/minmain/chapter2_9.html
 MIGRAD

This is the best minimizer for nearly all functions. It is a variable-metric method with inexact line search, a stable metric updating scheme, and checks for positive-definiteness. It will run faster if you [SET STRategy]SET STRATEGY 0 and will be more reliable if you [SET STRategy]SET STRATEGY 2 (although the latter option may not help much). Its main weakness is that it depends heavily on knowledge of the first derivatives, and fails miserably if they are very inaccurate. If first derivatives are a problem, they can be calculated analytically inside FCN (see elsewhere in this writeup) or if this is not feasible, the user can try to improve the accuracy of Minuit's numerical approximation by adjusting values using the [SET EPSmachine]SET EPS and/or [SET STRategy]SET STRATEGY commands (see Floating Point Precision and [SET STRategy]SET STRATEGY). 

http://wwwasdoc.web.cern.ch/wwwasdoc/WWW/minuit/minmain/chapter2_8.html
 SET STRategy<level>

Sets the strategy to be used in calculating first and second derivatives and in certain minimization methods. In general, low values of <level> mean fewer function calls and high values mean more reliable minimization. Currently allowed values are 0, 1 (default), and 2. 

http://wwwasdoc.web.cern.ch/wwwasdoc/WWW/minuit/minmain/chapter2_8.html
 MIGrad[maxcalls] [tolerance]

Causes minimization of the function by the method of Migrad, the most efficient and complete single method, recommended for general functions (see also MINImize). The minimization produces as a by-product the error matrix of the parameters, which is usually reliable unless warning messages are produced. The optional argument [maxcalls] specifies the (approximate) maximum number of function calls after which the calculation will be stopped even if it has not yet converged. The optional argument [tolerance] specifies required tolerance on the function value at the minimum. The default tolerance is 0.1, and the minimization will stop when the estimated vertical distance to the minimum (EDM) is less than 0.001*[tolerance]*UP (see [SET ERRordef]SET ERR). 

http://wwwasdoc.web.cern.ch/wwwasdoc/minuit/node18.html
 IMProve 	[maxcalls]
If a previous minimization has converged, and the current values of the parameters therefore correspond to a local minimum of the function, this command requests a search for additional distinct local minima. The optional argument [maxcalls] specifies the (approximate) maximum number of function calls after which the calculation will be stopped. 

http://wwwasdoc.web.cern.ch/wwwasdoc/WWW/minuit/minmain/chapter2_9.html
 MINIMIZE

This is equivalent to [MIGrad]MIGRAD, except that if [MIGrad]MIGRAD fails, it reverts to [SIMplex]SIMPLEX and then calls [MIGrad]MIGRAD again. This is what the old [MIGrad]MIGRAD command used to do, but it was removed from the [MIGrad]MIGRAD command so that users would have a choice, and because it is seldom of any use to call [SIMplex]SIMPLEX when [MIGrad]MIGRAD has failed (there are of course exceptions). 
   */

  const Double_t tol=  fTolerance;

  Double_t arglist[2];
  Int_t ierr = -999;

  //======================================================
  //higher stra, more difficult convergence
  //for each guided convergence, start with highest strategy with initial value from last convergent results, which were obtained with highest possible strategy

  Int_t stra = -999;
  Int_t s0 = 1;
  Int_t s1 = 1;
  if(klong){
    s0 = 2;
    s1 = 0;
  }

  for(stra = s0; stra>=s1; stra--){
    printf("-----------> MN2D::DoMigrad tol %e loop in stra %d startnow ...\n", tol, stra); 
    
    //always start from fMinuit=0x0, delete previous
    delete fMinuit;
    fMinuit = new TMinuit(Npar());
    
    fMinuit->SetFCN(MinFCNChisquare);
    fMinuit->SetPrintLevel(0);
    
    fMinuit->SetErrorDef(0.5);//ML 
    
    fStopwatch->Start(kTRUE);
    
    for(Int_t i=0;i<Npar();i++){
      const Double_t start=fPars[i];
      const Double_t ss=1e-2;
      const Double_t low= fLows[i];
      const Double_t hi= fHighs[i];
      
      fMinuit->DefineParameter(i, Form("p%d",i), start, ss, low,hi);
      if(fIfixed[i]){
        fMinuit->FixParameter(i);
      }
    }
    if(fMinuit->GetNumFreePars()!=Nfree()){
      printf("MN2D::DoMigrad fMinuit->GetNumFreePars()!=Nfree() %d %d\n", fMinuit->GetNumFreePars(), Nfree());exit(1);  
    }
    
    arglist[0]=stra;
    fMinuit->mnexcm("SET STRATEGY", arglist, 1, ierr);
    
    arglist[0] = MaxIter();
    arglist[1] = tol;

    if(klong){
      //MINIMIZE is much better than MIGRAD -- no, SIMPLEX introduces funny bias!
      fMinuit->mnexcm("MINIMIZE", arglist, 2, ierr);
    }
    else{
      fMinuit->mnexcm("MIGRAD", arglist, 2, ierr);
    }

    printf("-----------> MN2D::DoMigrad nbeta %03d tol %e finish ierr %d STRATEGY %d\n", nbeta, tol, ierr, stra);

    if(ierr==0){
      break;
    }
  }

  fStopwatch->Stop(); 
  fkErr = ierr;
  fStopwatch->Print("m"); //"m": milli sec; "u": micro sec

  return 0;
}

void MN2D::BetaFit(Int_t & nbeta, const Bool_t kHESSE)
{
  //======================================================================================
  //               minimization in each step (nbeta) of teh guided convergence
  //               called by ::Fit
  //               calling further sub function
  //               ::DoMigrad
  //======================================================================================

  //do migrad with global fBeta
  //fPars only modified after ::DoMigrad

  //update crossing using new fPars
  SetCrossing();

  //fix low mom p k fraction 0 and remained fixed
  FixOverflow();

  if(fMinuit){
    printf("MN2D::BetaFit exit fMinuit already exist!!! It should only be initialized in DoMigrad!!\n"); exit(1);
  }
  
  Status(Form("-------> MN2D::BetaFit nbeta %d startnow", nbeta));
  
  //every time start from fPars which is not updated since entering BetaFit()
  DoMigrad(nbeta, kHESSE);
  
  //now modify fPars
  //get par right after MIGRAD before HESSE modifies them
  Double_t dummy=-999;
  for(Int_t ipar=0;ipar<Npar();ipar++){
    fMinuit->GetParameter(ipar,fPars[ipar], dummy);
  }

  //update fStatChi2 and fPenalty with fPars
  GetFCNChisquare(fPars);
  Status(Form("-------> MN2D::BetaFit nbeta %d afterUpdatefPars", nbeta));

  //save intermediate results and output to file!

  if(kHESSE){
    //important!!
    //==> with the penalty defined as delta/StatErr, the total FCN has probability interpretation and the TMinuit err has correct meaning
    Int_t hesseflag = -999;
    Double_t hessearg[1]={MaxIter()};
    fStopwatch->Start(kTRUE);

    //to get rid of the influence by weighting!!
    //trigger true statistics
    //can not use true statistics because it is already resampled!!
    //--> must use so that if the sample size is double, stat. error will decrease
    fkUseTrueYield = kTRUE;

    fMinuit->mnexcm("HESSE", hessearg ,1, hesseflag);
    printf("-------> MN2D::BetaFit nbeta %d hesse done. flag %d\n", nbeta, hesseflag);

    //deactivate
    fkUseTrueYield = kFALSE;

    fStopwatch->Stop();
    fStopwatch->Print("m"); //"m": milli sec; "u": micro sec   
  }
  else{
    printf("-------> MN2D::BetaFit nbeta %d no hesse -- not needed in intermediate fits\n", nbeta);
  }

  for(Int_t ipar=0;ipar<Npar();ipar++){
    Double_t tmperr = -999, dummy = -999;
    fMinuit->GetParameter(ipar, dummy, tmperr);
    fLowErrs[ipar]=tmperr;
    fHighErrs[ipar]=tmperr;
  }
  TMatrixDSym *covmat = BaseUtils::GetCov(fMinuit, fIfixed, fLowErrs);

  SetPhysics(covmat, kHESSE);

  //only with kHESSE=1 the errors are correct!
  Save(Form("outmn2d_post_kHESSE%d_fkErr%d_iBeta%03d", kHESSE, fkErr, nbeta));

  nbeta++;

  delete covmat;
  delete fMinuit;
  //important: after delete, fMinuit can still have value. If later there is another delete, e.g. in Clear(), then it crashes!!
  fMinuit = 0x0;
}

void MN2D::Fit(const Int_t startnbeta)
{
  //===============================================
  //       Main coherent fit entry here!
  //       calling subfitting function
  //       ::BetaFit 
  //       which saves the output every time
  //       IMPORTANT: reduced chisuqare is important, normally if larger than 1.x, then the fitting fails even it reports convergence. Though this is rare -- ~1% chance.
  //===============================================

  fkUseTrueYield = kFALSE;

  Save("outmn2d_iBetaPre");

  Status(Form("---> MN2D::Fit() startnow"));

  Int_t nbeta = startnbeta;

  //fix all dEdx parameter beyond basic ones
  if(dEdxGausN::IsSignalShift()){
    fNfix += dEdxGausN::SetIfixExtraBB(fIfixed, kTRUE);
  }

  //======================================================
  //fit without regularization

  fBeta = 0;
  BetaFit(nbeta, kFALSE);

  //======================================================
  //enter regularization, fNPolint!=0 : has penality

  if(fNPolint){
    if(fPenalty<EPSILON){
      printf("MN2D::Fit exit penalty 0 !! %e when fNPolint %d\n", fPenalty, fNPolint);exit(1);
    }
    
    //in case beta!=0 (i.e. with regularization), release all limits
    for(Int_t ipar=0; ipar<Npar(); ipar++){
      fLows[ipar] = 0;
      fHighs[ipar]= 0;
    }

    //release all dEdx parameter beyond basic ones and fit once more
    //fNfix += dEdxGausN::SetIfixExtraBB(fIfixed, kFALSE);

    //======================================================
    //starting guided convergence, determine the lambda_0 (here is fBeta)
    const Double_t defactor = 2;
    Double_t kpenal = fPenalty;
    while(fStatChi2 > 0.5*kpenal){
      kpenal*=2;
      Status(Form("---> MN2D::Fit moderate penalty, increasing kpenal to %e", kpenal));
    }
    //very large penalty at first
    //start from chi2 ~ beta*penalty ==> beta<1
    fBeta = TMath::Power(defactor, (Int_t)(TMath::Log(fStatChi2/kpenal)/TMath::Log(defactor)));
    Status((Char_t*)"---> MN2D::Fit very large penalty");

    //guided convergence, increasing fBeta
    while( fBeta<=1){
      if(nbeta>10)
        break;

      BetaFit(nbeta, kFALSE);

      Status(Form("---> MN2D::Fit InLoopAfterBetaFit nbeta %d", nbeta-1));
      
      fBeta *= defactor;
    }

    //======================================================
    //finishing...

    fBeta = 1;
    //standard mean, res fixed from now
    //00x, 01x
    nbeta+=100;
    FinalFit(nbeta);
  }
  //===============================================

  //just to store pars
  Double_t standardpars[Npar()];
  memcpy(standardpars, fPars, sizeof(Double_t)*Npar());


  if(dEdxGausN::IsSignalShift()){
    nbeta+=100;
    //-0.025 even more far away --> too far away, data can not make such a long journey, go back to -0.02
    
    //for 10d+e pp corr=-0.02: data -> 0, MC -> -0.01 -- both far from initial value
    DoShift(standardpars, nbeta, 2,  fCpIni);

    
    nbeta+=100;
    //use fit value from previous fit
    DoShift(fPars, nbeta, 2, 0);

  }
}

void MN2D::FinalFit(Int_t &nbeta)
{
  if(fkErr){
    BetaFit(nbeta, kFALSE);
    Status(Form("---> MN2D::FinalFit refit due to previous failure beta %d ierr = %d output as iBeta%03d", nbeta-1, fkErr, nbeta-1));
  }

  //fix if not done yet
  //including also extra BB, 
  fNfix += dEdxGausN::SetIfixAllBB(fIfixed, kTRUE);
  //fixing res leads to convergence difficulty in shift proton -> No, floating makes it worse
  fNfix += SetIfixAllRes(fIfixed, kTRUE);

  //only call HESSE when mean res are fixed!
  nbeta+=10;
  BetaFit(nbeta, kTRUE);
  Status(Form("---> MN2D::FinalFit last Beta nbeta %d ierr = %d successful, copied as iBeta%03d", nbeta-1, fkErr, nbeta-1));

}

void MN2D::DoShift(const Double_t standardpars[], Int_t &beta, const Int_t itarget, const Double_t corr)
{
  //start from standard parameters

  //shifts are already 0 in standardpars
  memcpy(fPars, standardpars, sizeof(Double_t)*Npar());

  fPars[dEdxGausN::NparBasicBB() + itarget]= corr; 

  //_________________________________________________________________________________________
  //_________________________________________________________________________________________

  //still fixed, no HESSE
  if(TMath::Abs(corr)>EPSILON){
    BetaFit(beta, kFALSE);
    Status(Form("---> MN2D::FixMeanFit beta %d ierr = %d output as iBeta%03d", beta-1, fkErr, beta-1));
  }
  else{
    beta++;
  }

  //_________________________________________________________________________________________
  //_________________________________________________________________________________________

  //float, no HESSE
  fIfixed[dEdxGausN::NparBasicBB() + itarget]  = kFALSE; fNfix--;
  BetaFit(beta, kFALSE);
  Status(Form("---> MN2D::FixMeanFit beta %d ierr = %d output as iBeta%03d", beta-1, fkErr, beta-1));

  //_________________________________________________________________________________________
  //_________________________________________________________________________________________

  //fixed with HESSE
  //fix all extra
  FinalFit(beta);
}


Double_t MN2D::ScanShift(const Double_t standardpars[], Int_t &beta, const Int_t itarget, const Double_t p1, const Double_t p2, const Double_t dp)
{
  memcpy(fPars, standardpars, sizeof(Double_t)*Npar());

  for(Int_t ii=0; ii< dEdxGausN::NparSignalShift(); ii++){
    fPars[dEdxGausN::NparBasicBB()+ii] = 0;
  }

  TGraph * grScan = new TGraph;
  grScan->SetName(Form("grScan%d",itarget));

  //___________________________________________________________________________
  //___________________________________________________________________________

  Double_t minfcn = 1e10;
  Double_t bestScan = -999;
  Int_t bestBeta = -999;

  const Double_t pstart = dp>0? TMath::Min(p1, p2) : TMath::Max(p1, p2);
  const Double_t pstop  = dp>0? TMath::Max(p1, p2) : TMath::Min(p1, p2);

  for(Double_t pscan = pstart; (dp>0? pscan<=pstop:  pscan>=pstop) ; pscan+=dp){
    fPars[dEdxGausN::NparBasicBB() + itarget] = pscan;

    //all still fixed                                                                                                                                                                                                                                                                                         
    BetaFit(beta, kFALSE);
    Status(Form("---> MN2D::FixMeanFit beta %d itarget %d pscan %e ierr = %d output as iBeta%03d", beta-1, itarget, pscan, fkErr, beta-1));

    const Int_t ng = grScan->GetN();
    grScan->SetPoint(ng, pscan, myFCN());

    if(myFCN() < minfcn){
      minfcn=myFCN();
      bestScan = pscan;
      bestBeta = beta-1;
    }
  }

  printf("---> MN2D::ScanShift bestBeta %d bestScan %e minfcn %e\n", bestBeta, bestScan, minfcn);

  TFile * fout = new TFile(Form("outScan%d_bestBeta%03d_bestScan%.3f_minfcn%.0f.root", itarget, bestBeta, bestScan, minfcn), "recreate");
  grScan->Write();
  fout->Save();
  fout->Close();
  delete fout;

  delete grScan;

  return bestScan;
}

void MN2D::ScaleAll(const Double_t standardpars[], Int_t &beta, const Double_t scale)
{
  //start from standard parameters
  memcpy(fPars, standardpars, sizeof(Double_t)*Npar());

  //all mean and rms and shift fixed
  fPars[0] *= 1+scale;

  BetaFit(beta, kTRUE);
  Status(Form("---> MN2D::FixMeanFit beta %d ierr = %d output as iBeta%03d", beta-1, fkErr, beta-1));
}

